knitr::opts_chunk$set(echo = TRUE,
warning = FALSE,
message = FALSE)
library(tidyverse) # Default.
library(here) # To create file paths to data within different folders.
library(sf) # To read in spatial data and make simple features.
library(raster) # To work with raster data.
library(fasterize) # To make simple features and geometry.
NOTE: Working with raster data and images on the R Bren Server 3.6 will produce chunky, pixelated raster images on the .Rmd, but will produce clean, smooth raster images when knitted to an html.
knitr::include_graphics(here("img",
"landsat.png"))
landsat_file <- here("data",
"Landsat7.tif")
ls_1 <- raster(landsat_file)
ls_1
## class : RasterLayer
## band : 1 (of 5 bands)
## dimensions : 1758, 3701, 6506358 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : -59564.57, 51465.43, -404675.9, -351935.9 (xmin, xmax, ymin, ymax)
## crs : +proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0
## source : /home/csimms/esm_244_adv_data_analysis/labs/lab_6/data/Landsat7.tif
## names : Landsat7
## values : 0, 255 (min, max)
plot(ls_1)
ls_2 <- raster(landsat_file,
band = 2)
ls_3 <- raster(landsat_file,
band = 3)
ls_4 <- raster(landsat_file,
band = 4)
ls_stack <- raster::stack(landsat_file) # To read all the layers in at once.
ls_stack
## class : RasterStack
## dimensions : 1758, 3701, 6506358, 5 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : -59564.57, 51465.43, -404675.9, -351935.9 (xmin, xmax, ymin, ymax)
## crs : +proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0
## names : Landsat7.1, Landsat7.2, Landsat7.3, Landsat7.4, Landsat7.5
## min values : 0, 0, 0, 0, 0
## max values : 255, 255, 255, 255, 255
ls_1 <- raster::aggregate(ls_1,
fact = 3,
fun = mean)
ls_2 <- raster::aggregate(ls_2,
fact = 3,
fun = mean)
ls_3 <- raster::aggregate(ls_3,
fact = 3,
fun = mean)
ls_4 <- raster::aggregate(ls_4,
fact = 3,
fun = mean)
plot(ls_1,
col = hcl.colors(n = 100,
palette = 'Blues 2')) # To base a color scheme on a single palette.
plot(ls_2,
col = hcl.colors(n = 100,
palette = 'Greens 2'))
plot(ls_3,
col = hcl.colors(n = 100,
palette = 'Reds 2'))
plot(ls_4,
col = hcl.colors(n = 100,
palette = 'Reds 2'))
sbc_rast <- raster(here("data",
"county.tif"))
plot(sbc_rast)
plot(ls_3)
# ERROR Code in Lab Playthrough:
# mask(ls_3, sbc_rast) %>% plot()
# ls_3 <- mask(ls_3, sbc_rast)
# ls_4 <- mask(ls_4, sbc_rast)
vec1 <- 1:5
vec1
## [1] 1 2 3 4 5
vec1*2
## [1] 2 4 6 8 10
vec1^2
## [1] 1 4 9 16 25
ls_3
## class : RasterLayer
## dimensions : 586, 1234, 723124 (nrow, ncol, ncell)
## resolution : 90, 90 (x, y)
## extent : -59564.57, 51495.43, -404675.9, -351935.9 (xmin, xmax, ymin, ymax)
## crs : +proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0
## source : memory
## names : Landsat7
## values : 8.444444, 255 (min, max)
ls_3*2
## class : RasterLayer
## dimensions : 586, 1234, 723124 (nrow, ncol, ncell)
## resolution : 90, 90 (x, y)
## extent : -59564.57, 51495.43, -404675.9, -351935.9 (xmin, xmax, ymin, ymax)
## crs : +proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0
## source : memory
## names : Landsat7
## values : 16.88889, 510 (min, max)
log(ls_3)
## class : RasterLayer
## dimensions : 586, 1234, 723124 (nrow, ncol, ncell)
## resolution : 90, 90 (x, y)
## extent : -59564.57, 51495.43, -404675.9, -351935.9 (xmin, xmax, ymin, ymax)
## crs : +proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0
## source : memory
## names : layer
## values : 2.133509, 5.541264 (min, max)
plot(ls_3); plot(log(ls_3)) # To compare.
vec2 <- 6:10
vec1+vec2
## [1] 7 9 11 13 15
ls_3+ls_4
## class : RasterLayer
## dimensions : 586, 1234, 723124 (nrow, ncol, ncell)
## resolution : 90, 90 (x, y)
## extent : -59564.57, 51495.43, -404675.9, -351935.9 (xmin, xmax, ymin, ymax)
## crs : +proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0
## source : memory
## names : layer
## values : 35.33333, 498.7778 (min, max)
ls_stack <- stack(ls_1,
ls_2,
ls_3,
ls_4)
ls_mean <- raster::calc(ls_stack,
fun = mean,
na.rm = FALSE)
plot(ls_mean)
knitr::include_graphics(here("img",
"spectrum.png"))
knitr::include_graphics(here("img",
"ir_photo.jpg"))
\[NDVI = \frac{NIR - Red}{NIR + Red}\]
ndvi <- (ls_4 - ls_3) / (ls_4 + ls_3)
plot(ndvi,
col = hcl.colors(100,
palette = 'Grays'))
is_forest <- function(x,
thresh = 0.3) {
y <- ifelse(x >= thresh,
1,
NA)
return(y)
}
forest <- calc(ndvi,
fun = is_forest)
plot(forest,
col = 'green')
ndvi_df <- raster::rasterToPoints(ndvi) %>%
as.data.frame()
forest_df <- raster::rasterToPoints(forest) %>%
as.data.frame()
ggplot(data = ndvi_df,
aes(x = x,
y = y,
fill = layer)) +
geom_raster() +
geom_raster(data = forest_df,
fill = "green") +
coord_sf(expand = 0) +
scale_fill_gradient(low = "black",
high = "white") +
theme_void() +
theme(panel.background = element_rect(fill = "slateblue4"))